Phase resetting of collective rhythm in ensembles of oscillators 
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Phase resetting curves characterize the way a system with a collective periodic behavior responds 
to perturbations. We consider globally coupled ensembles of Sakaguchi-Kuramoto oscillators, and 
use the Ott-Antonsen theory of ensemble evolution to derive the analytical phase resetting equations. 
We show the final phase reset value to be composed of two parts: an immediate phase reset directly 
caused by the perturbation, and the dynamical phase reset resulting from the relaxation of the 
perturbed system back to its dynamical equilibrium. Analytical, semi-analytical and numerical 
approximations of the final phase resetting curve are constructed. We support our findings with 
extensive numerical evidence involving identical and non-identical oscillators. The validity of our 
theory is discussed in the context of large ensembles approximating the thermodynamic limit. 



I. INTRODUCTION 

The behavior of many biological systems displays 
rhythmic patterns that can be globally described through 
a periodic phase variable [l|. Examples are found on a 
variety of time-scales and levels of complexity: pulsat- 
ing neurons, circadian clocks in living beings, seasonal 
dynamics etc. [iHsf. The response of a periodic system 
to stimuli can be characterized by measuring the phase 
shift occurring as a result of a stimulus, i.e. by quanti- 
fying the difference between the collective phases of the 
perturbed and the original system. The phase shift's de- 
pendence on the phase value at which the perturbation 
occurred is termed phase resetting curve (or phase re- 
sponse curve, PRC), and is an inherent characteristic of 
any oscillatory system [J, Q . Study and applications of 
PRCs arc currently receiving a growing attention through 
their theoretical and experimental aspects [2, Q • 

Periodic biological systems are often composed of many 
interacting units, whose collective behavior emerges from 
the dynamical properties of the single units [2| . The col- 
lective (macroscopic) phase response of an empirical os- 
cillator is generated from stimuli that act on the level 
of individual (microscopic) units, and composed of their 
individual phase responses. 

In the context of neural oscillations various methods 
have been developed for the empirical measurement of 
PRCs, both in vivo and in vitro. The methods depend 
on the type/size of the neural ensemble, and are typically 
involving electric stimulation of neurons that induce a 
measurable phase shift in their spiking behavior 0, H 0] ■ 
Recently, the epileptiform activity in rats was character- 
ized through experimental PRCs, obtained by stimulat- 
ing thalamus and cortex of epileptic animals during spike 
and wave seizures Q . Investigation of PRCs is relevant 
for understanding the interaction properties of the neu- 
ral networks, such as their stability [7[, or synchroniza- 
tion and clustering [9|. The assumption of weak inter- 
neuron coupling was experimentally tested via infinites- 
imal PRCs [lOl- A new technique of PRC construction 
from data, based on weighted spike-triggered averaging 
was recently proposed lll[. While mostly studied in the 
domain of neurons, PRCs are also explored in other natu- 



ral rhythmic systems, that can be modeled as ensembles 
of many oscillatory elements. Empirical and theoreti- 
cal examination was thus conducted in coupled circadian 
clocks of insects [l^l , and the periodically driven saline 
oscillator [l3| . 

Due to its global periodic behavior, dynamics of neu- 
rons/oscillators and the resulting PRCs can be easily 
studied analytically or modeled computationally 0, M, 
Il4| . On the level of few neurons, extensive analytical re- 
sults including bifurcation diagrams are available for var- 
ious neuron models and types of interaction [3, [3, [iJ, [lal . 
On the other hand, global cooperative behavior in large 
neuron ensembles is typically examined numerically by 
modeling single elements as simple phase-oscillators. In 
this context, Kuramoto-type oscillators [J| are usually 
employed for their elegant and useful synchronization 
properties [l6| , that are well understood for various topo- 
logically and structurally different networks [13, [3 ■ ^ 
recent study of the phase response of stochastic oscil- 
lators revealed the interplay between microscopic and 
macroscopic collective phase sensitivity, for ensembles of 
globally coupled [1^ , and network coupled elements [201 ■ 
Phase resetting in globally coupled ensembles of oscilla- 
tors was investigated in relation to the symmetry prop- 
erties of the coupling function [2l| , emphasizing the pe- 
culiarity of non-odd coupling functional forms. PRCs 
were also constructed for models involving pulse-coupled 
neurons, providing insights into their interaction de- 
tails and collective dynamics at small |22| and large 
scales [2j| . Phase response behavior of complex oscilla- 
tor networks, besides giving insights into their interaction 
patterns [3, [2J] , can also assist in revealing the details of 
their network connectivity |25| . 

While the PRCs of a single (or few) oscillators are very 
well understood [j, |y, 0, [T5| , the results for ensembles of 
oscillators arc still rather scarce and often limited to nu- 
merical investigations only [3, [20, l21| . Namely, the diffi- 
culty here is that there are two contributions to the phase 
resetting. Apart from the phase shift immediately caused 
by the perturbation, the relaxation of the ensemble into 
a stationary state may induce an additional (potentially 
significant) phase shift, since any finite stimulus in gen- 
eral perturbs the system out of its equilibrium. 



In this paper wc consider a globally coupled ensemble 
of Sakaguchi-Kuramoto oscillators [26[, interacting via 
a non-odd coupling function. The system is perturbed 
from its dynamical equilibrium by applying the kick (of 
controllable form and strength) to the individual oscil- 
lators. We employ the Watanabe-Strogatz ansatz [24I 
in the recently proposed Ott-Antonsen (OA) formula- 
tion [28J, and calculate analytically the effect of the re- 
laxation of system's collective phase. As we show, an 
additional phase reset is induced, and its properties can 
be obtained from the ensemble's state measured immedi- 
ately after the perturbation. We confirm our theoretical 
results through a series of numerical simulations, consid- 
ering both identical and non-identical oscillators. The ex- 
tent of our theory is illustrated by considering the phase- 
resetting in an ensemble of Stuart-Landau oscillators. 

In opposition to [i^,l20|, here we consider deterministic 
oscillators and perturbations of arbitrary strength (not 
infinitesimal) , thus allowing the exploration of the phase 
reset's time evolution. 

This paper is organized as follows: in Section |TT] we 
formulate our model of kicked oscillator ensemble, define 
the concept of phase resetting and describe our imple- 
mentation of perturbation. In Section IIIII we show some 
analytical results regarding calculation of the immediate 
PRC. In Section [IVI we derive our main analytical result 
based on OA theory, for a general case of non-identical 
oscillators. We test our theory in Section |V] by exposing 
our numerical results for various oscillator ensembles and 
types. We discuss our results and conclude in SectionlVll 



II. FORMULATION OF OUR MODEL 

In this Section we construct our system of globally cou- 
pled oscillators and define the concept of phase resetting 
curve for the ensemble. We describe the implementation 
of the perturbation (kick) in response to which the PRC 
is measured. 



A. Ensemble and its description 

We consider an ensemble of N oscillators characterized 
by their natural frequencies ujk, whose dynamical states 
are defined by their phase values (pk G [0, 2tt). Oscillators 
are globally coupled using Sakaguchi-Kuramoto scheme 
with the interaction strength e > and phase shift /? G 



follows: 
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sm{(pj-ipk + l3)+q{ipk;ak) 5{t~T) . 

(1) 

The oscillator ipk undergoes a kick at time t = T, whose 
properties are given through the function q, which in gen- 
eral depends on the phase ipk and a set of parameters 
summarized in Sk- We take q to have a simple form as 



q{ip; a) ~ q{ip; A,a) ^ A sm{ip + a) 



(2) 



where A is the kicking strength, and a a phase shift pa- 
rameter. The collective dynamical state of the ensemble 
is quantified through the complex mean field (complex 
order parameter) Z defined as: 



1 

N 



k=l.N 



(e"^) = Re' 



where (•) stays for the average over ensemble. We call 
its absolute value R = \Z\ the collective radius, and its 
argument $ = arg Z the collective phase. The ensemble's 
PRC is defined with respect to the collective phase $. 

We now construct an alternative formulation of our 
system's equations by transforming the phase tp into a 
unitary complex variable a: 

a = e"^ . 

Thus, the time evolution of our kicked system given by 
Eq. Ilj now reads: 

&k = i(Tk Im iiujk + ee^'^Zal + r{ak]ak) 5{t - T)\ (3) 

where a* is the complex conjugate of a. The functional 
expression for the kick r becomes: 

r{a;d) ^r{a;A,\) = A\a (4) 

where A = e*". The complex mean field Z reads: 

k=l.N 

We term this approach unitary complex representation, in 
contrast to previously exposed phase representation. It 
transforms the trigonometric expressions in our equations 
into the algebraic ones, which will be of importance in our 
analytical studies that follow. Throughout this paper we 
will be interchangeably using both representations. 

Two types of phase-oscillator ensembles in respect to 
the frequency distributions are considered: 

{i) identical oscillators ujk = w, whose final dynamical 
state is the full synchronization pk = ^ (regardless of TV, 
e > and /3 G (— ^, ■!)). We will examine examples of 
small ensembles N < 00, and large ensembles approxi- 
mating thermodynamic limit TV — >■ 00. 

(ii) non-identical oscillators with Lorentzian frequency 
distribution g, characterized by a mean Co and a width 
7> 0: 



5(w) = - 



7 



n (uj — cD)2 4- 7^ 



(5) 



After relaxation this ensemble is partially synchronized 
(depending on (3 and 7). Here we consider only the ther- 
modynanuc limit iV — > cxd. The case of identical oscilla- 
tors is obtained for 7—^0+. 



B. General definition of PRC for the ensemble 

Consider an ensemble of TV > 1 oscillators governed by 
Eq. (H)) or ([3]) whose collective dynamical state is given 
by the complex mean field Z . We compare two realiza- 
tions of the ensemble, termed "original" and "kicked" 
system, with time evolutions described by Zit) and Z(t) 
respectively (all quantities related to the kicked system 
are denoted with bar). We assume the kicked system 
to be created at time t = T = in its initial state 
Z(0) = (a-(O)) = Zo = ((To) as the consequence of the 
kick acting on the original system, which at time t = 
is in its stationary state ^(0) = Zq: 



kick 



-^Zn 



or 



($o,i?o 



kick 



-^($o,i?o) 



After the kick, two systems independently continue their 
time evolutions - original system starting from Zq and 
kicked system starting from Zq. 

The phase resetting curve is defined as the difference 
between the collective phases of two systems: 

Z 

A = $ - $ = arg — 

/j 

as function of the collective phase value (of the original 
system) at which the kick occurred A ~ A($o)- We 
differ between two possible PRCs, depending on the time 
of phase reset measurement: 

• We call immediate phase resetting curve (pPRC) 
Aq = Ao($o) the phase shift value observed im- 
mediately after the resetting, i.e. as the immediate 
consequence of the kicking: 



Ao = $0 - *o = arg —- = arg 

^0 






(6) 



(we chose the abbreviation pPRC for "prompt 
PRC", since the abbreviation IPRC is in use for 
"infinitesimal PRC" O). 

• In contrast to pPRC, we term final phase resetting 
curve (fPRC) Aqo = Aoo(<I'o) the eventual phase 
reset value, measured when the kicked system re- 
laxes to its final stationary state: 

A^ = lim ($ - $) (t) = lim arg ^ . 

t-S-OO t-^QC Z(t) 

What we show in this paper, is that the curves Aq and 
Aoo arc in general different. We illustrate this on a sim- 
ple example. Consider iV = 10 identical oscillators whose 
collective dynamics is given by Eq. ((T]). They define 
the original system, synchronized {Rq = 1) at the phase 
value $0 where the kick occurs. The kick acts asymmet- 
rically: 5 out of 10 oscillators are kicked with a uniform 
strength A. The kicked system is created in its initial 
state Zq = Rqc^'^" with ^o < 1- while the original sys- 
tem continues to rotate with radius R{t) = 1. In Fig.[T]we 
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FIG. 1: (Color online) Ensemble of A*' = 10 identical os- 
cillators cijfe = Lj = 1 undergoing a kick at collective state 
Zo = e'*° (large dot). 5 oscillators are kicked with strength 
A = —0.25 (other oscillators are left unperturbed), for e = 
0.15 and fi = ^. Kicked system starts from its initial state 

Zo = -Roe'*" (large square). Time-evolutions of both complex 
mean fields are shown for comparison: original system Z(€) 
(dark/blue line) and kicked system Z(f) (light/red line). The 
values of $o — 'I'o = Aq, and A(t) are shown. 



show the evolution of the original system's complex mean 
field (blue) , in comparison with the one for the kicked sys- 
tem (red). The value of pPRC Aq = $o — ^o is clearly 
visible and can be calculated from Eq. ^ (we discuss this 
in the next Section). Aqo is measured once the kicked 
system reaches the synchronized state {R{t) — > 1), i.e. it 
is the limit value of A(t). As we show in what follows, 
the system's transient desynchronization induces differ- 
ent rotation speeds of the two collective phases which 
ultimately leads to a discrepancy between Aq and Aqo- 



C. Implementation of the kick 

We now construct the transformation that maps the 
phase of an oscillator from its pre-kick to its post-kick 
value: 

kick 

V'o > <^o 

in accordance with Eq. ^ and ^ . For iV = 1 and T = 
the Eq. ([T|) reduces to: 

1^ = w + j4sin(<y9 + a) 5{t) . 

We solve this equation by assuming that at the time t = 
the frequency term is negligible. This gives: 



Q(<Po) - Qifo) = 1 where Q 



dip 



^sin((/3 + a) 



which yields the relation 

tan^i±^=e^tan^i±^, 
2 2 

from which the final transformation is obtained: 



ifiQ = 2 arctan ( e tan ■ 



(7) 



We use the Eq. ([7]) for shifting the phase value of each 
oscillator in the ensemble as the consequence of the kick. 
Also, this expression automatically gives the PRC of an 
isolated oscillator A"^ in response to the kick Eq. ([2]): 



A-^ 



(^(0+) -1^(0 ) = (^0 -y'o • 



Since the properties of the time evolution of an isolated 
oscillator are not altered by the perturbation (no relax- 
ation back to stationary state) , its pPRC is always iden- 
tical to its fPRC. 

Alternatively, the kick transformation Eq. ([7|) can be 
formulated in terms of the unitary complex representa- 
tion. Proceeding equivalently as above, we obtain: 



where h = %-tt 

e^ + l 



Co 



tanh ■ 



Acto 



A - A26CT0 ' 



Furthermore, if we introduce 



the complex parameter tj ~ bX, the transformation above 
can be further simplified as: 



o-Q 



c^o -^7 



1 - ?/(To 

The PRC for a single oscillator A'^ now reads: 



(8) 



O-Q 



A — arg — = arg 



1 - iwoY 



1 



Wo 



The kick transformation given by Eq. (I7|) or (l8|) can be 
used together with Eq. ^ to analytically calculate Zq 
(and hence pPRC) for deterministic kicks. 



III. ANALYTICAL CALCULATION OF Zq 

In this Section we discuss the calculation of Zq in the 
thermodynamic limit for the cases when the kicking pa- 
rameters are picked from a given probability distribution. 
For a general ensemble with frequency distribution g, the 
complex mean field reads: 

Z~ dojdf g{oj)U{Lp,ijj)e'^ = / duida g{ijj)W{(7,u!)(T 

where U{ip^uj) is the fraction of oscillators having fre- 
quency w and phase ip (respectively, W{a, to) for w and 
a). After the kick, the complex mean field Zq = ((To) is: 

Zo = / dAda duj dipo g(w)p(w; A, a)U{ipo, w)e'^«('^°) 

= / dAda duj da^ g{uj)p{uj; A, a)W{cro, aj)a'o(cro) 

(9) 



where p stays for a given distribution of kicking param- 
eters A and a, that may depend on w. Thus, Zq can in 
principle always be calculated if the distributions U (or 
W) and p are known, and the pPRC can be obtained 
through Eq. ([6]) (although, the integral Eq. ([9]) can of- 
ten be solved only numerically). As our first numerical 
example, we study in Section |V] an ensemble of identical 
oscillators with kicking strengths A picked from a uni- 
form distribution A G [— Ao,Ao]. We compare the sim- 
ulation results with the 'analytical' curves obtained by 
calculating integrals Eq. © numerically. 

Alternatively, one can expand the expression Eq. (|S]) 
in power series to obtain: 



Za = ((o-o - '7*)Er=o(Wo) 

>o E,T=o(wo)") - Of Er=o(wo)' 



valid for |?7tTo| < 1 that is always verified in our case. We 
employ the expression above to directly calculate Zq for 
other two examples numerically studied in Section |Vl 

As the first of them, we consider identical oscillators 
with a unique kicking strength b = tanh ^ , and phase 
shifts a picked from the following distribution: 



p{a) = -- (1 + 2S'cos(a - a)) 

2n 



(10) 



Zo = (1 - 62)(e'*" + SbXe'^'^") - SbX 



where S < ^ and a S (— -j, ^) are real parameters. This 
distribution is characterized by (e*") = 5*6*" = SX and 
(e"""} = for all |m| > 1, i.e. only its first harmonic is 
non-zero. This means that (r/) = SbX, while (t?™) = for 
all \m\ > 1. Putting this in the series expression above, 
we obtain a closed simple formula for Zq: 

(11) 

As our final example, we examine non-identical oscilla- 
tors with Lorentzian frequency distribution Eq. ([S]). We 
assume the parameters b and A to be constant (but since 
the oscillators are not synchronized prior to the kick, this 
situation is not trivial). We have: 

where the quantities (Zq)^ = (ctq") are defined as in j28j . 
As shown there, the collective parameters describing the 
emergent dynamics of this system lie on OA manifold (we 
discuss this in detail in Section ITV| . defined by: 



(Zo)„ = (Zo)" for ah n>l. 

We assume that the departure from OA manifold due 
to the kicking is not significant, i.e. the identity above 
holds. This assumption gives: 



Zq = {--v*) 
V 



E 

n=0 



V-{Zor-- = ^^ (12) 

V 1 - vZo 



which is the OA approximation for Zq. Note that this is 
a direct generalization of the phase jump formula for a 
single oscillator Eq. ([5]). 



IV. TIME EVOLUTION OF THE PHASE RESET on the knowledge of Zq. 



In this Section we employ the Ott-Antonsen theory and 
estimate the additional phase shift due to the evolution 
after the kick. We construct a formula for fPRC Aqo for 
a general Sakaguchi-Kuramoto ensemble. 



A. An illustrative example 

We begin by illustrating the concepts of pPRC and 
fPRC further. Consider the simplest nontrivial ensem- 
ble consisting of two oscillators, with dynamics given by 
Eq. ([T]) or (j3]). We apply two different versions of the 
kicking in the system's synchronized state: 

(a) both oscillators kicked with the equal strengths Ai = 
A2 = 0.18 

(b) one oscillator kicked with Ai = 0.2, the other left 
unperturbed {A2 ~ 0) 

For both cases we numerically compute Aq and Aqq • Re- 
sults are reported in Fig.[2J while in the case (a) we have 
a perfect agreement between Aq and Aqo, in the case (b) 
fPRC Aoo significantly differs from the pPRC Aq. The 
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FIG. 2: (Color online) PRC for an ensemble of two identical 
oscillators having uj = 1, with two different kicking schemes, 
(a): Ai ^ A2 ^ 0.18, (b) Ai = 0.2, A2 = 0. Both cases: 
/3 = ^ and e — 0.1. Legend indicates various cases/curves. 



difference between the curves Aq and Aqo in case (b) is 
not uniform, but it depends on the collective phase value. 
In the case (a) both oscillators receive equal kicks and 
hence the same immediate phase shift. This leaves the 
system in the synchronized state i?o = ^0 = 1, generat- 
ing no transient relaxation: fPRC is given by the pPRC. 
In contrast to this, in the case (b) non-uniform kicking 
desynchronizes the system Rq < Rq, which generates an 
additional phase shift during the ensemble's relaxation 
(cf. Fig.[T|). In the reminder of this Section we provide a 
theoretical framework for approximating the fPRC based 



B. Ott-Antonsen theory for non-identical 
oscillators 



The Sakaguchi-Kuramoto ensemble of globally cou- 
pled identical oscillators can be exactly described by the 
Watanabe-Strogatz (WS) theory [23|. However, the evo- 
lution of mean field in this theory is described through 
the constants of motion, that are determined from the 
initial conditions by rather complicated expressions. Re- 
cently, Ott and Antonsen found a simple ansatz that 
leads to a closed set of equations for the mean field [231 ■ 
As shown in [29|, this ansatz corresponds to a special 
choice of constants of motion in the WS theory. Such a 
dynamical state (defining the OA manifold) is eventually 
reached in an ensemble of non-identical oscillators with 
Lorentzian frequency distribution. In this paper we ap- 
ply the OA theory by assuming that the deviations from 
OA manifold due to the kick are small (cf. derivation of 
the Eq. (|12p V We obtain simple (although approximate) 
analytical expressions, whose precision is shown through 
numerical results in Section IVl 

We thus consider the general case of an ensemble with 
collective dynamics given by Eq. ([1]) or ([3]), and frequen- 
cies UJ picked from the Lorentzian distribution Eq. ([S]). In 
the OA approximation, the evolution of the mean field Z 
(in thermodynamic limit) is given by [28|: 



Z = iojZ -jZ+-Z{e 



^/3 e-»/3|^|2N 



(13) 



After substituting Z = _Re'* this equation is decomposed 
into two real ODEs: 

■ e cos B „, „„, „ 
R = -—J1r{i^r^)-jR, 

■ esin/3. ^„, 
$ = 0J+ ^ (1 + fl^) • 



The only stable fixed point for the first equation is 

The system achieves this fi- 



Rf = Jl- 



27 



< 1. 



e cos [3 

nal radius value only if ecos/3 > 27, and otherwise re- 
mains fully dcsynchronized with R{t) — ?> (note that 
full synchronization Rf = 1 occurs only for 7 = 0). 
For R = Rf the collective phase rotates with frequency 
il = u! + e sin/3 — 7tan/3. Both equations above can be 
easily integrated, yielding: 



R{t) 
$(t) 

tan [3 
2 



Rfil 



Ri-R, 



2. -(e cos 0-27) 



•) 



-1/2 



= ^Q + ujt + e s'm (3 t - J tan /3 t+ 



where Rq = ^(0) and $0 = ^(0). Consider again the 
original and kicked ensemble. The original system is 
at time i = in stationary state Zq, with radius Rq 



(equals to final stationary radius Rq ~ Rf) and phase 
evolution $(i) = $o + ^t- The kicked system starts at 
time t = from the initial state Zq = /Joe**", with the 
time evolution given by the equation above. We con- 
sider the phase difference between two systems defined 
as A(i) = ($ - $)(t) which reads: 



A(i) = (do - $0 



i^Mlnd 



"lz^p-(ecos/3-27) t 



)+tan/3 1n(|^ 



(14) 

Coming back to the definitions of pPRC and fPRC from 
Section ini we differ between two situations: 

• at t = last two terms in Eq. p^ cancel out, 
leaving 

A(f = 0) = $0 - $0 = Ao 

which is the pPRC Ag, consistently with what ex- 
posed in Section im 

• as f — > cxD the middle term on RHS of Eq. (fH|) 
vanishes at the limit, leaving 



Aoo 



Ao + A 



R 



where we have defined 



tan/3 In 



i?o 

Rn 



tan/3 In 



Zo 
Zo 



(15) 



(note that Rq = Rf = Rf, as the final stationary 
i?- value depends only on the ensemble parameters) . 

The term A^ is what accounts for the additional phase 
shift due to system's relaxation back to stationary state. 
It depends on the nature of perturbation through -^ and 
the dynamics of kicked system through /3. For perturba- 
tions that leave Rq = Rq (e.g. equal perturbation applied 
to all oscillators as in Fig. [5] case (a)) and systems with 
/3 = 0, the fPRC Aoo is identical to pPRC Aq. On the 
other hand, if /3 7^ the rotation speed of the collective 
phase depends on the value of R. Hence, if the system's 
stationary radius is perturbed {Rq ^ Rq), an additional 
phase reset A/j is generated at the limit (cf. Fig. [T] and 
Fig- m case (b)), which depends on the kicking phase $o- 
If the details of the kick are known An can be cal- 
culated from Eq. (|15p and the analytical approximation 
for fPRC can be constructed. Alternatively, fPRC can 
be numerically approximated by summing the computed 
curves Aq and A^, even without knowing the parame- 
ters that specify the system and the kick. We stress again 
that Eq. p^ is approximate, as it is based on the OA 
ansatz; its precision is related to the deviations from the 
OA manifold. 



C. Ott-Antonsen theory for identical oscillators 

The equivalent results for the ensemble of identical 
oscillators can be easily recovered from what just ex- 
posed by taking 7 — )■ 0+. Since identical oscillators are 



always synchronized prior to the kicking, we also have 
Rq = Rf = Rf = 1, and Rq < 1- The phase difference 
equation Eq. ([Tl)) reduces to: 

A(i)= ($o-*o)+ 

^ln(l + i^e-"°^^*)+tan/3 lni?o, 

from which we obtain the formula for fPRC: 

Aoo = Aq + tan /? In | Zq | . 

The dynamics of the complex order parameter given by 
Eq. (jl3|) only weakly depends on 7, and the limit 7 ^■ 
0^ in this equation is not singular. However, the whole 
validity of the OA theory leading to Eq. p^ heavily 
depends on the spreading range of the frequencies 7, and 
the limit 7 — !■ 0"*" app ears to be singular in this respect 
(see discussion in [28ll29|). 

In the next Section we show the numerical results con- 
firming the theory just exposed. 



V. NUMERICAL RESULTS 

In this Section we expose our numerical results, con- 
firming our theoretical findings from previous Sections. 
For a given ensemble consisting of N oscillators, we start 
from a random initial phase for each oscillator ip G [0, 27r) 
and run the dynamics until the system reaches the sta- 
tionary state. Then, a collective phase value $0 is chosen 
at random (in the stationary regime) and the kick is ap- 
plied (in accordance with Eq. (O for each oscillator) . The 
phase reset is immediately measured, and the respective 
values of Ao($o) and Aij($o) are recorded. The kicked 
system is then relaxed into a stationary state, when the 
final phase reset Aoo($o) is measured with respect to 
the original system (at the same final time value). Re- 
peating this procedure for many values of $0 yields the 
numerical curves Aq, A.^ and Aoo. These are compared 
with the corresponding analytical curves obtained as pre- 
viously discussed. In the last paragraph we examine the 
phase- resetting in an ensemble of complex Stuart-Landau 
oscillators, showing the validity our results in a broader 
context. 



A. Small ensembles of identical oscillators 

We start with an ensemble of A^ = 10 identical oscil- 
lators with the kick implemented as follows: 3 oscilla- 
tors are kicked with a strength oi Ai, 4 are kicked with 
a strength of A2, while the remaining 3 are left unper- 
turbed. For this case Zq can be easily calculated using 
Eqs.® & dH), which gives Aq and A^, and hence the 
fPRC Aoo • In Fig. [3] we show the computed numerical 
curves Aq, A/j and Aoo (marked by different symbols), in 
comparison with the analytical fPRC Aoo (thick curve). 
The OA equation Eq. (jl3|) is an approximation valid at 
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FIG. 3: (Color online) Ensemble of A'' = 10 identical oscil- 
lators, 3 kicked with a strength Ai — 0.05, 4 kicked with a 
strength A2 — —0.06, and the remaining 3 unperturbed. Pa- 
rameters: Qfc = 0, uj = 1, e — 0.1, P = ^. Numerical and 
analytical curves are indicated in Legend. 



the thermodynamic limit. Nevertheless, in Fig. [3] analyt- 
ical fPRC approximates the numerical fPRC with a good 
precision, despite small ensemble size. Wc considered 
many different deterministic kick realizations (through 
parameters A and a) for TV = 10 identical oscillators, 
and in general obtained a very good agreement between 
the numerical and analytical fPRC. Note that the mag- 
nitude of Afl is comparable to that of Aq - the post-kick 
evolution of the phase reset is significant. 

We now pick the kicking parameters at random: A^- 
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FIG. 4: (Color online) Ensemble of A*' = 10 identical os- 
cillators with random kicks: strengths Ak picked from a 
Gaussian distribution with mean and standard deviation 
0.1, Qffc-values picked uniformly from [— 7r,7r]. Parameters: 
UJ = 1, £ = 0.1, j3 — ^. Numerical curves shown as indi- 
cated in Legend. 



values from a Gaussian distribution centered at 0, and 
afc-values uniformly from interval [— 7r,7r]. Numerically 
computed curves Aq, A/j and A 00 are reported in Fig. |4] 
(same symbols are used). Since the details of the kick are 
now unknown, we are unable to construct the analytical 
approximations. However, it follows from the formula 
Eq. (|15p that one can sum numerically obtained Aq and 
A/j and approximate Aoo- We show this in Fig. |3J nu- 
merical curves (Aq -|- A^j) and Aoo agree rather well, thus 
confirming the practical validity of the formula Eq. (|15p . 



B. Large ensembles of identical oscillators 

Wc now examine large ensembles approximating the 
thermodynamic limit through three examples discussed 
in Section Hill We start with N = 1000 identical oscilla- 
tors and the kick implemented by uniformly picking the 
kicking strengths A from the interval [— Ao,Ao], while 
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FIG. 5: (Color online) Ensemble oi N = 1000 identical os- 
cillators with kicking strengths A picked from A £ [—^0,^0] 
for Ao = 0.1, with aU q = 0, and lo = I, e = 0.1, /3 = ^. 
All the curves are shown as indicated in the Legend. Top: a 
single realization of A-values. Bottom: 10 values taken from 
each of 60 runs with different realizations of A- values. 



setting a — ior all oscillators. We compute the nu- 
merical curves Ag, Aji and Aqo, and construct the cor- 
responding 'analytical' curves from Eq. ([9]) (quotes in- 
dicate that the integral Eq. ^ was solved numerically). 
Both numerical and analytical results are shown in Fig. [5] 
(same symbols as previously). The top figure reports 
results of a run with a single realization of A-values, 
while in the bottom figure we show sampled data ob- 
tained from many different realization of A-values. For 
completeness, we now report all three analytical curves 
(differently dashed thick lines). There is a very good 
agreement between all numerical and analytical curves, 
even for a single realization of A- values, due to the large 
ensemble size. For the sampled curves, note that the 
theoretical fPRC almost perfectly approximates the av- 
eraged value of numerical fPRC. This further confirms 
the quality of the OA approximation. 

We now set all the kicking strengths to constant A = 
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FIG. 6: (Color online) Ensemble of iV = 1000 identical os- 
cillators. All kicking strengths are fixed to Ao = 0.1, while 
a-values are picked from [— 7r,7r] via distribution Eq. (|10|l for 
5 = 0.15 and a = TT. oj = 1, e = 0.1, l3 = ^. AH the 
curves are shown as indicated in the Legend. Top: a single 
realization of a-values. Bottom: 10 PRC-values taken from 
each of 60 runs with different realizations of a-values. 



Aq, while picking the a-values from the interval [— 7r,7r] 
via the distribution Eq. (jlOp . We use Eq. ([TT|) to test our 
numerical findings. Note that in opposition to the pre- 
vious example, now we are having an entirely analytical 
solution for the whole PRC-problem. The numerical and 
analytical results are reported in Fig. |6l where we again 
show a case of a single a-values realization (top) , and the 
sampled data from many different realization of a-values 
(bottom). We have a very good agreement between all 
obtained curves, particularly in the sampled case. The 
agreement is however not perfect, and it can be used to 
estimate the precision of OA approximation. 



C. Large ensembles of non-identical oscillators 

We now examine an ensemble consisting of A^ = 1000 
non-identical oscillators with a Lorentzian frequency dis- 
tribution Eq. ([5]). All kicking strengths are set to A = Aq, 
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FIG. 7: (Color online) Ensemble of A = 1000 non-identical 
oscillators with Lorentzian frequency distribution Eq. ([5]) with 
LJ = 10. All kicking strengths y4o = 0.1 and a = 0. Other 
parameters: e = 1, /3 = ^. Top: 7 = 0.2, bottom: 7 = 
0.3. All the curves are shown as indicated in the Legend. 



Frequencies cj^ are taken to be ujk — 7tan( 



with a = 0. Since the ensemble is not fuUy synchronized 
prior to the kicking, phase jumps of osciUators are in 
general different, depending on their phases value at the 
time of kick. We now use Eq. p^ to obtain the an- 
alytical curves, assuming the kick does not remove the 
system significantly from the OA manifold. The results 
are shown in Fig. [7] for two different values of Lorentzian 
width 7. We show only the cases of single kicking re- 
alizations (since they are always noisy). Again, a very 
good agreement between all the numerical and analytical 
curves is obtained, even in the case of 7 = 0.3 (bottom) 
where we observe a discontinuity in the pPRC and fPRC. 
This confirms the assumption of small departure from OA 
manifold due to kicking, used to obtain Eq. p^ . 
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D. Large ensembles of Stuart-Landau oscillators 

In the final part of this Section, we wish to illustrate 
the extent of the theory developed in this work. We 
consider an ensemble of A^ = 1000 non-identical Stuart- 
Landau oscillators |j|, described by the complex ampli- 
tudes w = re^^ . The system's equation which reads: 



w 



(e + lo; - i\w\'^)w + ee'^Z - A'e-*"'(5(i) 



can be decomposed into: 



r = ^(1 — r 



ip: 



^ eR cos{^ - If + /3) 

-A'cos{ip + a') S{t) 

w + ef sin($ -ip + (3) + ^ sin((y9 + a') 6(1) 



where Z = i?e** = (w) is the mean field, through which 
the ensemble is coupled. For ^ 3> 1 the system's behavior 
resembles phase-oscillators since r w 1 after transients 
for all oscillators. This makes the non-kicked part of the 
phase equation above reduce to the Eq. ([T|). The kick 
is implemented in analogy with the Section |lTl but now 
acts according to: 

fg = Iroe'vo -A'e'"'\ 

ipo = arg(roe"^" -A'e'"') 

In equivalence with the previous paragraph, we consider 
an ensemble with a single realization of the Lorentzian 
frequency distribution, and a uniform kick A' ~ A'q, a' ~ 
0. In Fig.|S]we show the numerical curves, along with the 
analytical ones which are obtained from the Eq. ([T2|) . but 
using A — 0.022 instead of A ^ 0.1. Namely, the kick 
of a strength A'q = 0.1 described by the equations above 
produces different shift in the collective phase, roughly 
five times smaller than a kick of strength ^ = 0.1 given by 
Eq. ([7]) (since it affects both r and (p). Nevertheless, there 
is a good agreement between numerical and analytical 
curves, despite having considered a completely different 
(complex, i.e. two-dimensional) class of oscillators. This 
suggests that the range of application of both OA theory 
and our phase-resetting theory extends beyond simple 
phase-oscillators. 



FIG. 8: (Color online) Ensemble of iV = 1000 Stuart- 
Landau oscillators with Lorentzian frequency distribution 
(i = 10, 7 = 0.2 (cf. Fig. [71). All kicking strengths 
A'o = 0.1, a' = 0. Other parameters: e = 1, fi ^ ^, ^ = 10. 
All the curves are shown as indicated in the Legend (analyti- 
cal ones are obtained using Eq. (|12|l for A = 0.022 - see text 
for details). 



VI. CONCLUSIONS 

We constructed a theory of the phase resetting of the 
collective phase for ensembles of Sakaguchi-Kuramoto os- 
cillators interacting globally via a non-odd coupling func- 
tion. By employing the Ott-Antonsen theory we showed 
that the final phase resetting curve can be recovered from 
the collective system's state Zq, measured immediately 
after the perturbation. We confirmed our theory on a 
series of numerical examples involving small and large 
ensembles of identical and non-identical oscillators, us- 
ing various kicking patterns and strengths. In particular, 
we showed an example where the entire PRC-problem 
was solved analytically (Eq. pT|)L Finally, by investigat- 
ing PRC of an ensemble of Stuart-Landau oscillators, we 
showed that our theory is valid beyond phase-oscillator 
models. 

Our main result captured by formula Eq. (jlSp , provides 
three methods of obtaining the final PRC Aqo : 

• numerically, from (a long) run (dark/blue symbols 
in Figs, nil), 

• analytically, by calculating the curve Aqo using 
Eqs.dS]) & (US]) (full thick curve in Figs. [H]), 



• semi- analytically, by applying the Eq. (jisp to the 
numerically obtained curves Aq and An (squares 
in Fig.lH). 

Our numerical results show a good agreement for all three 
methods on all examined examples. 

Our findings suggest that the system's relaxation af- 
ter the perturbation can significantly alter its final PRC. 
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In agreement with [2l|, we found that the initial and 
the final phase resets coincide for odd interaction func- 
tions {(3 ~ 0) or perturbations that leave the collective 
state unchanged {Rq = Rq). However, once the collec- 
tive radius R is perturbed, different rotation speeds of the 
collective phases of the original and kicked systems (pro- 
portional to sin/3 ^ 0) induce an additional non-trivial 
phase reset at the limit. 

It is worth noticing the case of \(3\ = ^: here our the- 
ory no longer holds, since the mean field is zero Z = 0, 
and therefore the collective phase is not well defined. 
This renders impossible to generally define the concept 
of phase-resetting. However, in this case the kick tran- 
siently induces Z ^ contrasting the disordered behav- 
ior for Z = - this situation can perhaps be used to 
model the dynamical behavior underlying the appearance 
of Event Related Potentials, which are of g reat current 
interest in neuroscience and psychology [30|. 

Our theory is valid for arbitrary Sakaguchi-Kuramoto- 
type ensembles of identical oscillators (e.g., a nonlinear 
coupling |3l| can be easily incorporated), and ensembles 
of non-identical oscillators with Lorentzian frequency dis- 
tribution. Its limitations are in principle given by the 
range of validity of OA theory [2^, [2^ |31|; nevertheless, 
our numerical results indicate that our analytical findings 
apply even for conditions that are formally not included 
in the OA theory, such as small ensembles (far from ther- 
modynamic limit). 



An important extension of our results regards the gen- 
eralization on complex networks. A general network con- 
nection topology will generate a more complicated post- 
kick time evolution of the collective phase, which is not 
covered by the OA theory. Understanding evolution of 
the phase reset on complex networks is of high relevance 
for problems related to the reconstruction of the network 
structure from the experimental data. By investigating 
the properties of the empirical PRCs, one might be able 
to recover the connectivity patterns of the underlying 
network Ii,[2|. 

Another relevant question revolves around the stochas- 
tic periodic systems, and different methods of defining 
their PRCs [S^l ■ A complete theory of stochastic PRCs 
will be of great help in constructing the PRC from peri- 
odic (but noisy) experimental data. Finally, since a PRC 
can be generally defined for any rhythmic system, our 
results can be extended to other types/models of oscilla- 
tors, such as different models of neural cells. 
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